For an individual, the ability to cope with challenges is fundamental to surviving in a changing environment. The corticosterone stress response has a major impact on the ability of an organism to react to the challenges of rapid environmental change. Previous research has shown that when faced with a prolonged cold snap, adult female tree swallows increase their baseline and stress-induced corticosterone levels, and that individual variation in corticosterone levels predicts reproductive success during challenges. Tree swallows are advancing reproduction with warming springs and are more likely to encounter cold snaps during breeding, which can reduce reproductive success. Here, using data from a long-term study of free-living tree swallows collected from 2013-2025, we show that nestling tree swallows show an immediate increase in corticosterone levels when ambient temperatures fall below a threshold of 18.5°C, which has previously been identified as the threshold at which insect availability begins to decline. As a result of earlier reproduction, nestlings are more likely to experience cold snaps during development, making it especially important to understand their effects. Ongoing work is also testing whether male and female nestlings display dimorphism in their corticosterone responses to further our understanding of how nestlings are responding to thermal environmental shifts.
Hypothesis 1: Temperature will show a negative relationship with baseline and stress-induced corticosterone, with there being an increase in these measures when both average daily temperatures and average 3-day temperatures are lower; Temperature will have no relationship with post-dex corticosterone levels.
Hypothesis 2: Males and females will show base differences in their corticosterone regulation, with females showing lower corticosterone across all measures, based on previous studies done in adult birds (baseline corticosterone, stress-induced corticosterone, and post-dex corticosterone).
Hypothesis 3:Females will showcase a stronger corticosterone response (across all measures) in response to lower temperatures, consistent with prior work in our population showing that females may more sensitive to some stressors.
In this section, I am running GAMM models, as based on prior checking, the effect of temperature on corticosterone is non-linear, and a smooth term fits the data better than a linear term. Average daily temperature is used as the predictor, as this seemed to have a slightly stronger effect than average 3 day temperature, and all corticosterone measures have been log-transformed to conform to a normal distribution better.
bleed1_weather_gamm4 <- gamm4( # run gamm model
log_bleed1_final_cort ~ s(avgC_capture_day) + (bleed1_latency_sec),
random = ~(1 | nest_key),
data = weather_df
)
summary(bleed1_weather_gamm4$gam) # view results of gam
##
## Family: gaussian
## Link function: identity
##
## Formula:
## log_bleed1_final_cort ~ s(avgC_capture_day) + (bleed1_latency_sec)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.993948 0.114410 8.688 <2e-16 ***
## bleed1_latency_sec 0.001967 0.001082 1.818 0.0695 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(avgC_capture_day) 3.703 3.703 23.39 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.206
## lmer.REML = 1709 Scale est. = 0.73143 n = 601
summary(bleed1_weather_gamm4$mer) # view random effects
## Linear mixed model fit by REML ['lmerMod']
##
## REML criterion at convergence: 1709
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.8056 -0.6282 -0.0074 0.5488 3.5072
##
## Random effects:
## Groups Name Variance Std.Dev.
## nest_key (Intercept) 0.3886 0.6233
## Xr s(avgC_capture_day) 1.2794 1.1311
## Residual 0.7314 0.8552
## Number of obs: 601, groups: nest_key, 164; Xr, 8
##
## Fixed effects:
## Estimate Std. Error t value
## X(Intercept) 0.993948 0.113681 8.743
## Xbleed1_latency_sec 0.001967 0.001082 1.818
## Xs(avgC_capture_day)Fx1 -0.556697 0.323104 -1.723
##
## Correlation of Fixed Effects:
## X(Int) Xbl1__
## Xbld1_ltnc_ -0.841
## Xs(vgC__)F1 -0.029 0.032
plot(bleed1_weather_gamm4$gam, pages = 1) # view smooth term
k.check(bleed1_weather_gamm4$gam) # check the k to see if the model is "wiggly" enough
## k' edf k-index p-value
## s(avgC_capture_day) 9 3.703261 0.8236044 0
par(mfrow = c(2, 2))
gam.check(bleed1_weather_gamm4$gam) # check residuals to make sure model is normal
##
## 'gamm' based fit - care required with interpretation.
## Checks based on working residuals may be misleading.
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
##
## k' edf k-index p-value
## s(avgC_capture_day) 9.0 3.7 0.82 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
base_cort <- draw(bleed1_weather_gamm4$gam, term = "avgC_capture_day")$data %>%
dplyr::select(avgC_capture_day = avgC_capture_day, fit = .estimate, se = .se, cil = .lower_ci, ciu = .upper_ci)
base_cort_plot <- ggplot() +
geom_abline(slope = 0, intercept = 0, linetype = "61", color = "gray75", size = 0.3) +
geom_line(data = base_cort, aes(x = avgC_capture_day, y = fit), color = base_color, linetype = "solid", size = 0.3) +
geom_ribbon(data = base_cort, aes(x = avgC_capture_day, ymin = cil, ymax = ciu), fill = base_color, alpha = 0.4) +
labs(x = NULL, y = "Predicted Partial Baseline Corticosterone") +
theme_classic() +
scale_x_continuous(breaks = seq(0, 30, 10)) +
xlim(c(13, 30)) +
theme(panel.grid = element_blank(), axis.text = element_text(size = 12),axis.title = element_text(size = 14, family = "Helvetica"), plot.margin = margin(15, 15, 15, 15))
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.
base_cort_plot
bleed2_weather_gamm4 <- gamm4( # run gamm model
log_bleed2_final_cort ~ s(avgC_capture_day) + (log_bleed1_final_cort),
random = ~(1 | nest_key),
data = weather_df
)
summary(bleed2_weather_gamm4$gam) # view results of gam
##
## Family: gaussian
## Link function: identity
##
## Formula:
## log_bleed2_final_cort ~ s(avgC_capture_day) + (log_bleed1_final_cort)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.32155 0.07382 31.447 < 2e-16 ***
## log_bleed1_final_cort 0.22625 0.03584 6.313 5.54e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(avgC_capture_day) 3.544 3.544 6.734 0.000112 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.219
## lmer.REML = 1447.7 Scale est. = 0.48454 n = 572
summary(bleed2_weather_gamm4$mer) # view random effects
## Linear mixed model fit by REML ['lmerMod']
##
## REML criterion at convergence: 1447.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.1550 -0.5149 0.0832 0.5904 2.3502
##
## Random effects:
## Groups Name Variance Std.Dev.
## nest_key (Intercept) 0.4551 0.6746
## Xr s(avgC_capture_day) 1.2671 1.1256
## Residual 0.4845 0.6961
## Number of obs: 572, groups: nest_key, 158; Xr, 8
##
## Fixed effects:
## Estimate Std. Error t value
## X(Intercept) 2.32155 0.07431 31.242
## Xlog_bleed1_final_cort 0.22625 0.03606 6.275
## Xs(avgC_capture_day)Fx1 -0.45712 0.32667 -1.399
##
## Correlation of Fixed Effects:
## X(Int) Xl_1__
## Xlg_bld1_f_ -0.542
## Xs(vgC__)F1 -0.038 0.069
plot(bleed2_weather_gamm4$gam, pages = 1) # view smooth term
k.check(bleed2_weather_gamm4$gam) # check the k to see if the model is "wiggly" enough
## k' edf k-index p-value
## s(avgC_capture_day) 9 3.544486 0.6504017 0
par(mfrow = c(2, 2))
gam.check(bleed2_weather_gamm4$gam) # check residuals to make sure model is normal
##
## 'gamm' based fit - care required with interpretation.
## Checks based on working residuals may be misleading.
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
##
## k' edf k-index p-value
## s(avgC_capture_day) 9.00 3.54 0.65 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
stress_cort <- draw(bleed2_weather_gamm4$gam, term = "avgC_capture_day")$data %>%
dplyr::select(avgC_capture_day = avgC_capture_day, fit = .estimate, se = .se, cil = .lower_ci, ciu = .upper_ci)
stress_cort_plot <- ggplot() +
geom_abline(slope = 0, intercept = 0, linetype = "61", color = "gray75", size = 0.3) +
geom_line(data = stress_cort, aes(x = avgC_capture_day, y = fit), color = purple_color, linetype = "solid", size = 0.3) +
geom_ribbon(data = stress_cort, aes(x = avgC_capture_day, ymin = cil, ymax = ciu), fill = purple_color, alpha = 0.4) +
labs(x = NULL, y = "Predicted Partial Stress-Induced Corticosterone") +
theme_classic() +
scale_x_continuous(breaks = seq(0, 30, 10)) +
xlim(c(13, 30)) +
theme(panel.grid = element_blank(), axis.text = element_text(size = 12), axis.title = element_text(size = 14, family = "Helvetica"), plot.margin = margin(15, 15, 15, 15))
## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.
stress_cort_plot
bleed3_weather_gamm4 <- gamm4( # run gamm model
log_bleed3_final_cort ~ s(avgC_capture_day) + (log_bleed1_final_cort),
random = ~(1 | nest_key),
data = weather_df
)
summary(bleed3_weather_gamm4$gam) # view results of gam
##
## Family: gaussian
## Link function: identity
##
## Formula:
## log_bleed3_final_cort ~ s(avgC_capture_day) + (log_bleed1_final_cort)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.77662 0.06524 27.233 <2e-16 ***
## log_bleed1_final_cort -0.01010 0.03350 -0.301 0.763
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(avgC_capture_day) 2.852 2.852 4.569 0.00815 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.0417
## lmer.REML = 954.39 Scale est. = 0.32079 n = 454
summary(bleed3_weather_gamm4$mer) # view random effects
## Linear mixed model fit by REML ['lmerMod']
##
## REML criterion at convergence: 954.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.5213 -0.5110 -0.0397 0.5311 5.1130
##
## Random effects:
## Groups Name Variance Std.Dev.
## nest_key (Intercept) 0.2768 0.5261
## Xr s(avgC_capture_day) 0.3262 0.5711
## Residual 0.3208 0.5664
## Number of obs: 454, groups: nest_key, 121; Xr, 8
##
## Fixed effects:
## Estimate Std. Error t value
## X(Intercept) 1.77662 0.06500 27.333
## Xlog_bleed1_final_cort -0.01010 0.03337 -0.303
## Xs(avgC_capture_day)Fx1 0.24314 0.20741 1.172
##
## Correlation of Fixed Effects:
## X(Int) Xl_1__
## Xlg_bld1_f_ -0.516
## Xs(vgC__)F1 0.047 -0.096
plot(bleed3_weather_gamm4$gam, pages = 1) # view smooth term
k.check(bleed3_weather_gamm4$gam) # check the k to see if the model is "wiggly" enough
## k' edf k-index p-value
## s(avgC_capture_day) 9 2.851564 0.7074518 0
par(mfrow = c(2, 2))
gam.check(bleed3_weather_gamm4$gam) # check residuals to make sure model is normal
##
## 'gamm' based fit - care required with interpretation.
## Checks based on working residuals may be misleading.
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
##
## k' edf k-index p-value
## s(avgC_capture_day) 9.00 2.85 0.71 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
dex_cort <- draw(bleed3_weather_gamm4$gam, term = "avgC_capture_day")$data %>%
dplyr::select(avgC_capture_day = avgC_capture_day, fit = .estimate, se = .se, cil = .lower_ci, ciu = .upper_ci)
dex_cort_plot <- ggplot() +
geom_abline(slope = 0, intercept = 0, linetype = "61", color = "gray75", size = 0.3) +
geom_line(data = dex_cort, aes(x = avgC_capture_day, y = fit), color = dex_color, linetype = "solid", size = 0.3) +
geom_ribbon(data = dex_cort, aes(x = avgC_capture_day, ymin = cil, ymax = ciu), fill = dex_color, alpha = 0.4) +
labs(x = NULL, y = "Predicted Partial Post-dexamethasone Corticosterone") +
theme_classic() +
scale_x_continuous(breaks = seq(0, 30, 10)) +
xlim(c(13, 30)) +
theme(panel.grid = element_blank(), axis.text = element_text(size = 12), axis.title = element_text(size = 14, family = "Helvetica"), plot.margin = margin(15, 15, 15, 15),
)
## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.
dex_cort_plot
Combine all plots together
combined <- base_cort_plot | stress_cort_plot | dex_cort_plot # one row
final_figure1 <- combined + plot_annotation(
tag_levels = "A",
caption = "Average Daily Temperature (°C)"
) &
theme(
plot.tag = element_text(
size = 18,
family = "Helvetica",
face = "bold"
),
plot.tag.position = c(0.25, 0.98), # inside the panel (x,y in npc coords)
plot.caption = element_text(
hjust = 0.5,
size = 14,
family = "Helvetica"
)
)
final_figure1
ggsave("TablesFigures/Figure1.png",
plot = final_figure1, height = 5, width = 14)
Here, I am just using regular linear mixed models to determine if there are any sex differences.
bleed1_sex <- lmer(log_bleed1_final_cort ~ sex + bleed1_latency_sec + ( 1 | nest_key), data = sex_df)
summary(bleed1_sex)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log_bleed1_final_cort ~ sex + bleed1_latency_sec + (1 | nest_key)
## Data: sex_df
##
## REML criterion at convergence: 1757.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.4701 -0.5800 -0.0154 0.5371 3.4105
##
## Random effects:
## Groups Name Variance Std.Dev.
## nest_key (Intercept) 0.6510 0.8069
## Residual 0.7263 0.8522
## Number of obs: 601, groups: nest_key, 164
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 9.494e-01 1.274e-01 5.389e+02 7.454 3.63e-13 ***
## sexMale 3.989e-02 7.865e-02 5.153e+02 0.507 0.612
## bleed1_latency_sec 2.010e-03 1.103e-03 5.433e+02 1.822 0.069 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) sexMal
## sexMale -0.282
## bld1_ltncy_ -0.765 -0.003
simulationOutput = simulateResiduals(bleed1_sex, plot = F)
plot(simulationOutput)
testDispersion(simulationOutput)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 1.0048, p-value = 0.864
## alternative hypothesis: two.sided
bleed2_sex <- lmer(log_bleed2_final_cort ~ sex + log_bleed1_final_cort + ( 1 | nest_key), data = sex_df)
summary(bleed2_sex)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log_bleed2_final_cort ~ sex + log_bleed1_final_cort + (1 | nest_key)
## Data: sex_df
##
## REML criterion at convergence: 1459.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.1207 -0.5130 0.0789 0.5499 2.2702
##
## Random effects:
## Groups Name Variance Std.Dev.
## nest_key (Intercept) 0.5201 0.7212
## Residual 0.4863 0.6974
## Number of obs: 572, groups: nest_key, 158
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 2.25749 0.08172 295.17103 27.626 < 2e-16 ***
## sexMale 0.05062 0.06597 475.13834 0.767 0.443
## log_bleed1_final_cort 0.25653 0.03526 561.00350 7.276 1.16e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) sexMal
## sexMale -0.359
## lg_bld1_fn_ -0.464 -0.029
simulationOutput = simulateResiduals(bleed2_sex, plot = F)
plot(simulationOutput)
testDispersion(simulationOutput)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 1.0182, p-value = 0.824
## alternative hypothesis: two.sided
bleed3_sex <- lmer(log_bleed3_final_cort ~ sex + log_bleed1_final_cort + ( 1 | nest_key), data = sex_df)
summary(bleed3_sex)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log_bleed3_final_cort ~ sex + log_bleed1_final_cort + (1 | nest_key)
## Data: sex_df
##
## REML criterion at convergence: 957.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.5645 -0.5057 -0.0177 0.5134 5.0510
##
## Random effects:
## Groups Name Variance Std.Dev.
## nest_key (Intercept) 0.2921 0.5404
## Residual 0.3224 0.5678
## Number of obs: 454, groups: nest_key, 121
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 1.79434 0.07105 221.04399 25.256 <2e-16 ***
## sexMale -0.09935 0.06002 378.72536 -1.655 0.0987 .
## log_bleed1_final_cort 0.01622 0.03227 446.74111 0.503 0.6154
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) sexMal
## sexMale -0.394
## lg_bld1_fn_ -0.444 -0.022
simulationOutput = simulateResiduals(bleed3_sex, plot = F)
plot(simulationOutput)
testDispersion(simulationOutput)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.97126, p-value = 0.8
## alternative hypothesis: two.sided
Here, I will be using GAM agains to fit that smooth term effect of temp on cort, to stay consistent with the first set of models.
bleed1_int_gamm <- gamm4(
log_bleed1_final_cort ~ sex + s(avgC_capture_day, by = sex) + (bleed1_latency_sec),
random = ~(1 | nest_key),
data = sex_df)
summary(bleed1_int_gamm$gam)
##
## Family: gaussian
## Link function: identity
##
## Formula:
## log_bleed1_final_cort ~ sex + s(avgC_capture_day, by = sex) +
## (bleed1_latency_sec)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.966837 0.119479 8.092 3.35e-15 ***
## sexMale 0.050639 0.077758 0.651 0.5151
## bleed1_latency_sec 0.002023 0.001085 1.865 0.0627 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(avgC_capture_day):sexFemale 3.216 3.216 17.27 <2e-16 ***
## s(avgC_capture_day):sexMale 3.406 3.406 16.92 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.202
## lmer.REML = 1719.1 Scale est. = 0.73292 n = 601
summary(bleed1_int_gamm$mer)
## Linear mixed model fit by REML ['lmerMod']
##
## REML criterion at convergence: 1719.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.6269 -0.6366 0.0073 0.5611 3.4329
##
## Random effects:
## Groups Name Variance Std.Dev.
## nest_key (Intercept) 0.3902 0.6247
## Xr.0 s(avgC_capture_day):sexMale 1.4040 1.1849
## Xr s(avgC_capture_day):sexFemale 1.0624 1.0307
## Residual 0.7329 0.8561
## Number of obs: 601, groups: nest_key, 164; Xr.0, 8; Xr, 8
##
## Fixed effects:
## Estimate Std. Error t value
## X(Intercept) 0.966837 0.119426 8.096
## XsexMale 0.050639 0.077617 0.652
## Xbleed1_latency_sec 0.002023 0.001085 1.865
## Xs(avgC_capture_day):sexFemaleFx1 -0.376943 0.314524 -1.198
## Xs(avgC_capture_day):sexMaleFx1 -0.654523 0.351124 -1.864
##
## Correlation of Fixed Effects:
## X(Int) XsexMl Xbl1__ X(C__):F
## XsexMale -0.298
## Xbld1_ltnc_ -0.801 -0.004
## Xs(C__):FF1 -0.035 0.003 0.037
## Xs(C__):MF1 -0.013 0.002 0.015 0.165
plot(bleed1_int_gamm$gam, pages = 1) # view smooth term
k.check(bleed1_int_gamm$gam) # check the k to see if the model is "wiggly" enough
## k' edf k-index p-value
## s(avgC_capture_day):sexFemale 9 3.215974 0.819781 0
## s(avgC_capture_day):sexMale 9 3.405542 0.819781 0
par(mfrow = c(2, 2))
gam.check(bleed1_int_gamm$gam) # check residuals to make sure model is normal
##
## 'gamm' based fit - care required with interpretation.
## Checks based on working residuals may be misleading.
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
##
## k' edf k-index p-value
## s(avgC_capture_day):sexFemale 9.00 3.22 0.82 <2e-16 ***
## s(avgC_capture_day):sexMale 9.00 3.41 0.82 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
pred_base <- ggpredict(bleed1_int_gamm,
terms = c("avgC_capture_day [all]", "sex [Female, Male]"))
base_int <- ggplot(pred_base, aes(x = x, y = predicted, color = group, fill = group)) +
geom_abline(slope = 0, intercept = 0, linetype = "61", color = "gray75", size = 0.3) +
geom_line(size = 1.0) +
geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = 0.4, color = NA) +
labs(x = NULL, y = "Predicted Baseline Corticosterone") +
theme_classic() +
theme(
panel.grid = element_blank(),
axis.text = element_text(size = 12),
axis.title = element_text(size = 14, family = "Helvetica"), plot.margin = margin(15, 15, 15, 15),
legend.title = element_blank(),
legend.position = "top") +
scale_color_manual(values = c(Female = female_color, Male = male_color))+
scale_fill_manual(values = c(Female = female_color, Male = male_color))+
scale_x_continuous(breaks = seq(13, 30, 5))
base_int
bleed2_int_gamm <- gamm4(
log_bleed2_final_cort ~ sex + s(avgC_capture_day, by = sex) + (log_bleed1_final_cort),
random = ~(1 | nest_key),
data = sex_df)
summary(bleed2_int_gamm$gam)
##
## Family: gaussian
## Link function: identity
##
## Formula:
## log_bleed2_final_cort ~ sex + s(avgC_capture_day, by = sex) +
## (log_bleed1_final_cort)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.28829 0.07935 28.836 < 2e-16 ***
## sexMale 0.05430 0.06558 0.828 0.408
## log_bleed1_final_cort 0.23400 0.03558 6.576 1.1e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(avgC_capture_day):sexFemale 3.642 3.642 6.012 0.00044 ***
## s(avgC_capture_day):sexMale 2.279 2.279 4.951 0.00707 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.212
## lmer.REML = 1453.8 Scale est. = 0.48081 n = 572
summary(bleed2_int_gamm$mer)
## Linear mixed model fit by REML ['lmerMod']
##
## REML criterion at convergence: 1453.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.1192 -0.4925 0.0871 0.5539 2.4230
##
## Random effects:
## Groups Name Variance Std.Dev.
## nest_key (Intercept) 0.4614 0.6792
## Xr.0 s(avgC_capture_day):sexMale 0.2509 0.5009
## Xr s(avgC_capture_day):sexFemale 1.4752 1.2146
## Residual 0.4808 0.6934
## Number of obs: 572, groups: nest_key, 158; Xr.0, 8; Xr, 8
##
## Fixed effects:
## Estimate Std. Error t value
## X(Intercept) 2.28829 0.07985 28.657
## XsexMale 0.05430 0.06543 0.830
## Xlog_bleed1_final_cort 0.23400 0.03589 6.519
## Xs(avgC_capture_day):sexFemaleFx1 -0.48520 0.33575 -1.445
## Xs(avgC_capture_day):sexMaleFx1 -0.24713 0.19585 -1.262
##
## Correlation of Fixed Effects:
## X(Int) XsexMl Xl_1__ X(C__):F
## XsexMale -0.361
## Xlg_bld1_f_ -0.490 -0.033
## Xs(C__):FF1 -0.024 0.002 0.036
## Xs(C__):MF1 -0.053 -0.007 0.093 0.187
plot(bleed2_int_gamm$gam, pages = 1) # view smooth term
k.check(bleed2_int_gamm$gam) # check the k to see if the model is "wiggly" enough
## k' edf k-index p-value
## s(avgC_capture_day):sexFemale 9 3.642472 0.6378559 0
## s(avgC_capture_day):sexMale 9 2.279260 0.6378559 0
par(mfrow = c(2, 2))
gam.check(bleed2_int_gamm$gam) # check residuals to make sure model is normal
##
## 'gamm' based fit - care required with interpretation.
## Checks based on working residuals may be misleading.
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
##
## k' edf k-index p-value
## s(avgC_capture_day):sexFemale 9.00 3.64 0.64 <2e-16 ***
## s(avgC_capture_day):sexMale 9.00 2.28 0.64 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
pred_stress <- ggpredict(bleed2_int_gamm,
terms = c("avgC_capture_day [all]", "sex [Female, Male]"))
stress_int <- ggplot(pred_stress, aes(x = x, y = predicted, color = group, fill = group)) +
geom_abline(slope = 0, intercept = 0, linetype = "61", color = "gray75", size = 0.3) +
geom_line(size = 1.0) +
geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = 0.4, color = NA) +
labs(x = NULL, y = "Predicted Stress-Induced Corticosterone") +
theme_classic() +
theme(panel.grid = element_blank(),
axis.text = element_text(size = 12),
axis.title = element_text(size = 14, family = "Helvetica"), plot.margin = margin(15, 15, 15, 15),
legend.title = element_blank(),
legend.position = "top") +
scale_color_manual(values = c(Female = female_color, Male = male_color))+
scale_fill_manual(values = c(Female = female_color, Male = male_color))+
scale_x_continuous(breaks = seq(13, 30, 5))
stress_int
bleed3_int_gamm <- gamm4(
log_bleed3_final_cort ~ sex + s(avgC_capture_day, by = sex) + (log_bleed1_final_cort),
random = ~(1 | nest_key),
data = sex_df)
summary(bleed3_int_gamm$gam)
##
## Family: gaussian
## Link function: identity
##
## Formula:
## log_bleed3_final_cort ~ sex + s(avgC_capture_day, by = sex) +
## (log_bleed1_final_cort)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.812871 0.071565 25.332 <2e-16 ***
## sexMale -0.094864 0.060123 -1.578 0.115
## log_bleed1_final_cort -0.001366 0.033306 -0.041 0.967
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(avgC_capture_day):sexFemale 2.15 2.15 2.349 0.1174
## s(avgC_capture_day):sexMale 2.29 2.29 2.572 0.0979 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.0317
## lmer.REML = 963.09 Scale est. = 0.32267 n = 454
summary(bleed3_int_gamm$mer)
## Linear mixed model fit by REML ['lmerMod']
##
## REML criterion at convergence: 963.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.4835 -0.5153 -0.0258 0.5166 5.0105
##
## Random effects:
## Groups Name Variance Std.Dev.
## nest_key (Intercept) 0.2779 0.5272
## Xr.0 s(avgC_capture_day):sexMale 0.1564 0.3955
## Xr s(avgC_capture_day):sexFemale 0.1199 0.3463
## Residual 0.3227 0.5680
## Number of obs: 454, groups: nest_key, 121; Xr.0, 8; Xr, 8
##
## Fixed effects:
## Estimate Std. Error t value
## X(Intercept) 1.812871 0.070667 25.654
## XsexMale -0.094864 0.060046 -1.580
## Xlog_bleed1_final_cort -0.001366 0.033224 -0.041
## Xs(avgC_capture_day):sexFemaleFx1 0.121431 0.146830 0.827
## Xs(avgC_capture_day):sexMaleFx1 0.128565 0.165693 0.776
##
## Correlation of Fixed Effects:
## X(Int) XsexMl Xl_1__ X(C__):F
## XsexMale -0.391
## Xlg_bld1_f_ -0.460 -0.031
## Xs(C__):FF1 0.030 0.003 -0.086
## Xs(C__):MF1 0.043 0.028 -0.096 0.174
plot(bleed3_int_gamm$gam, pages = 1) # view smooth term
k.check(bleed3_int_gamm$gam) # check the k to see if the model is "wiggly" enough
## k' edf k-index p-value
## s(avgC_capture_day):sexFemale 9 2.150400 0.6986927 0
## s(avgC_capture_day):sexMale 9 2.290195 0.6986927 0
par(mfrow = c(2, 2))
gam.check(bleed3_int_gamm$gam) # check residuals to make sure model is normal
##
## 'gamm' based fit - care required with interpretation.
## Checks based on working residuals may be misleading.
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
##
## k' edf k-index p-value
## s(avgC_capture_day):sexFemale 9.00 2.15 0.7 <2e-16 ***
## s(avgC_capture_day):sexMale 9.00 2.29 0.7 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
pred_dex <- ggpredict(bleed3_int_gamm,
terms = c("avgC_capture_day [all]", "sex [Female, Male]"))
dex_int <- ggplot(pred_dex, aes(x = x, y = predicted, color = group, fill = group)) +
geom_abline(slope = 0, intercept = 0, linetype = "61", color = "gray75", size = 0.3) +
geom_line(size = 1.0) +
geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = 0.4, color = NA) +
labs(x = NULL, y = "Predicted Post-dexamethasone Corticosterone") +
theme_classic() +
theme(panel.grid = element_blank(),
axis.text = element_text(size = 12),
axis.title = element_text(size = 14, family = "Helvetica"), plot.margin = margin(15, 15, 15, 15),
legend.title = element_blank(),
legend.position = "top") +
scale_color_manual(values = c(Female = female_color, Male = male_color))+
scale_fill_manual(values = c(Female = female_color, Male = male_color))+
scale_x_continuous(breaks = seq(13, 30, 5))
dex_int
Combine all together
combined2 <- base_int | stress_int | dex_int # one row
final_figure2 <- combined2 + plot_annotation(
tag_levels = "A",
caption = "Average Daily Temperature (°C)"
) &
theme(
plot.tag = element_text(
size = 18,
family = "Helvetica",
face = "bold"
),
plot.tag.position = c(0.25, 0.85),
plot.caption = element_text(
hjust = 0.5,
size = 14,
family = "Helvetica"
))
final_figure2
ggsave("TablesFigures/Figure2.png",
plot = final_figure2, height = 5, width = 14)
Now, creating a figure to compare the smooths in the interaction models, based on this tutorial: this tutorial:https://fromthebottomoftheheap.net/2017/10/11/difference-splines-i/
Run function from above site:
smooth_diff <- function(model, newdata, f1, f2, var, alpha = 0.05,
unconditional = FALSE) {
fac <- as.character(newdata[[var]]) # added this to convert sex to a character to allow for comparision
xp <- predict(model, newdata = newdata, type = 'lpmatrix')
c1 <- grepl(f1, colnames(xp))
c2 <- grepl(f2, colnames(xp))
r1 <- fac == f1 # modified this
r2 <- fac == f2
## difference rows of xp for data from comparison
X <- xp[r1, ] - xp[r2, ]
## zero out cols of X related to splines for other lochs
X[, ! (c1 | c2)] <- 0
## zero out the parametric cols
X[, !grepl('^s\\(', colnames(xp))] <- 0
dif <- X %*% coef(model)
se <- sqrt(rowSums((X %*% vcov(model, unconditional = unconditional)) * X))
crit <- qt(alpha/2, df.residual(model), lower.tail = FALSE)
upr <- dif + (crit * se)
lwr <- dif - (crit * se)
data.frame(pair = paste(f1, f2, sep = '-'),
diff = dif,
se = se,
upper = upr,
lower = lwr)
}
Create dataframe to be fed into the function:
pdat_base <- expand.grid(
avgC_capture_day = seq(13.5, 28.5, length = 400),
sex = c("Female", "Male"),
bleed1_latency_sec = mean(sex_df$bleed1_latency_sec, na.rm = TRUE)
)
pdat_stress_dex <- expand.grid(
avgC_capture_day = seq(13.5, 28.5, length = 400),
sex = c("Female", "Male"),
log_bleed1_final_cort = mean(sex_df$log_bleed1_final_cort, na.rm = TRUE)
)
comp_base <- smooth_diff(bleed1_int_gamm$gam, pdat_base, 'Female', 'Male', 'sex')
comp_base$avgC_capture_day <- pdat_base$avgC_capture_day[pdat_base$sex == "Female"]
comp_stress <- smooth_diff(bleed2_int_gamm$gam, pdat_stress_dex, 'Female', 'Male', 'sex')
comp_stress$avgC_capture_day <- pdat_stress_dex$avgC_capture_day[pdat_stress_dex$sex == "Female"]
comp_dex <- smooth_diff(bleed3_int_gamm$gam, pdat_stress_dex, 'Female', 'Male', 'sex')
comp_dex$avgC_capture_day <- pdat_stress_dex$avgC_capture_day[pdat_stress_dex$sex == "Female"]
Create plots:
base_diff <- ggplot(comp_base, aes(x = avgC_capture_day, y = diff)) +
geom_ribbon(aes(ymin = lower, ymax = upper), alpha = 0.2, fill = purple_color) +
geom_line(size = 1, color = purple_color) +
geom_hline(yintercept = 0, linetype = "dashed", color = "gray") +
coord_cartesian(ylim = c(-1, 1)) +
labs(x = NULL, y = 'Estimated Difference in Baseline Corticosterone \n(Female - Male)') +
theme_classic() +
theme(panel.grid = element_blank(), axis.text = element_text(size = 12),axis.title = element_text(size = 14, family = "Helvetica"), plot.margin = margin(15,15,15,15))
base_diff
## Stress-induced corticosterone
stress_diff <- ggplot(comp_stress, aes(x = avgC_capture_day, y = diff)) +
geom_ribbon(aes(ymin = lower, ymax = upper), alpha = 0.4, fill = purple_color) +
geom_line(size = 1, color = purple_color) +
geom_hline(yintercept = 0, linetype = "dashed", color = "gray") +
coord_cartesian(ylim = c(-1, 1)) +
labs(x = NULL, y = 'Difference in Stress-Induced Corticosterone\n(Female - Male)') +
theme_classic() +
theme(panel.grid = element_blank(), axis.text = element_text(size = 12),axis.title = element_text(size = 14, family = "Helvetica"), plot.margin = margin(15,15,15,15))
stress_diff
dex_diff <- ggplot(comp_dex, aes(x = avgC_capture_day, y = diff)) +
geom_ribbon(aes(ymin = lower, ymax = upper), alpha = 0.4, fill = purple_color) +
geom_line(size = 1, color = purple_color) +
geom_hline(yintercept = 0, linetype = "dashed", color = "gray") +
coord_cartesian(ylim = c(-1, 1)) +
labs(x = NULL, y = 'Difference in Post-dexamethasone Corticosterone\n(Female - Male)') +
theme_classic() +
theme(panel.grid = element_blank(), axis.text = element_text(size = 12),axis.title = element_text(size = 14, family = "Helvetica"), plot.margin = margin(15,15,15,15))
dex_diff
Create one figure with all combined:
combined3 <- base_diff | stress_diff | dex_diff # one row
final_figure3 <- combined3 +
plot_annotation(
tag_levels = "A",
caption = "Average Daily Temperature (°C)"
) &
theme(
plot.tag = element_text(
size = 18,
family = "Helvetica",
face = "bold"
),
plot.tag.position = c(0.25, 0.98),
plot.caption = element_text(
hjust = 0.5,
size = 14,
family = "Helvetica"
)
)
final_figure3
ggsave("TablesFigures/Figure3.png",
plot = final_figure3, height = 5, width = 14)